#include <iostream>
#include <vector>
#include "../include/Eigen/Dense"
#include "../include/norm.hpp"
#include <math.h>

using namespace std;
using namespace Eigen;

#define pi 4 * atan(1)

class Function
{
public:
    Function() {}
    double operator()(double _x)
    {
        return exp(sin(_x));
    }
    double f_x(double _x)
    {
        return exp(sin(_x)) * cos(_x);
    }
    double laplace_f(double _x) // 负拉普拉斯算子
    {
        return exp(sin(_x)) * (sin(_x) - pow(cos(_x),2));
    }
};

int main()
{
    Function func;
    int dimention = 33;
    double h = 1.0 / 32;
    MatrixXd A = MatrixXd::Zero(33, 33);
    VectorXd vec_f = VectorXd::Zero(dimention);
    for (int i = 0; i < dimention; i++)
    {
        if ((i == 0) || (i == dimention - 1))
        {
            A(i, i) = 1;
            vec_f(i) = func(i * h);
        }
        else
        {
            A(i, i - 1) = -1.0;
            A(i, i + 1) = -1.0;
            A(i, i) = 2.0;
            vec_f(i) = func.laplace_f(i * h) * pow(h, 2);
        }
    }
    cout << A << endl;
    VectorXd sol = A.lu().solve(vec_f);
    Norm nor;
    VectorXd real = VectorXd::Zero(dimention);
    for (int i = 0; i < dimention; i++)
    {
        real(i) = func(i * h);
    }
    double res = nor.inftynorm(sol - real);
    cout << "x = [";
    for (int i = 0; i < dimention; i++)
    {
        cout << i * h << ",";
    }
    cout << "];" << endl;

    cout << "y = [";
    for (int i = 0; i < dimention; i++)
    {
        cout << real[i] << ",";
    }
    cout << "];" << endl;

    cout << "scatter(x,y, \".\");" << endl;
}